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1.  THE  PROBLEM 


The  means  and  mechanism  for  translating  lamina  level  properties  into  laminate  properties 
has  been  available  for  many  years  under  strongly  limited  conditions  loosely  referred  to  as 
those  of  a  "thin"  laminate.  Th.s  procedure  is  known  as  classical  lamination  theory,  and  it 
neglects  the  interface  conditions  between  individual  lamina  and  directly  sums,  by  algebraic 
formulas,  the  inplane  lamina  properties  to  obtain  the  integrated  inplane  laminate  properties. 
This  two-dimensional  (2-D)  procedure  is  widely  used  and  has  contributed  greatly  to  the 
effectiveness  of  fiber  composites.  Unfortunately,  the  2-D  procedure  provides  little  guidance  on 
how  to  proceed  in  the  much  more  complicated  three-dimensional  (3-D)  case.  However, 
Pagano  (1974)  has  provided  a  complete  and  general  3-D  lamination  procedure,  and  it  has 
been  implemented  by  Sun  and  Li  (1988).  The  2-D,  classical  lamination  procedure  is 
approximate  in  the  sense  that  it  implies  plane  stress  conditions.  It  is  the  intention  here  to 
develop  a  3-D  lamination  procedure  which  retains  the  essential  simplicity  of  the  2-D  procedure. 
Of  course,  it  will  not  be  possible  to  assume  plane  stress  conditions  in  the  3-D  case,  ano  some 
other  approximations  will  need  to  be  introduced  here  to  accomplish  ihe  objective. 

The  starting  point  in  the  materials  characterization  is  the  properties  state  at  the  uni¬ 
directional  lamina  level.  Taking  these  properties  to  be  characterized  by  a  state  of  transversely 
isotropic  symmetry  leads  to  the  properties  specification  through  five  independent  properties. 
These  five  independent  elastic  properties  can  be  specified  by  the  longitudinal  Young’s 
modulus,  E{,  meaning  the  modulus  in  the  fiber  direction,  the  longitudinal  Poisson’s  ratio,  v|t 
meaning  the  transverse  strain  response  when  the  composite  is  strained  in  the  longitudinal 
direction  to  determine  E)t  as  weli  as  the  transverse  Young's  modulus,  £t,  the  longitudinal  shear 
modulus.  Mj,  and  finally  the  transverse  shear  modulus,  pt.  These  five  lamina  level  properties 
can  be  differentiated  by  grouping  them  into  fiber-dominated  vs.  matrix-dominated  properties. 
The  two  grouping  are: 


Fiber  Dominated 


Matrix  Dominated 


For  the  fiber-dominated  properties,  the  fibers  act  as  the  direct  load  transfer  agent.  In  the 
matrix-dominated  properties,  the  fibers  effectively  act  as  "inclusions"  in  a  continuous  matrix 
phase.  This  distinction  is  fundamental,  and  reveals  itself  in  all  credible  micro-mechanics 
models  of  fiber  composites.  High  performance  fibers  translate  into  a  fiber-property  dominated 
composite.  Accordingly,  in  a  fiber-dominated  composite,  the  matrix-dominated  properties  are 
of  lesser  importance  than  are  the  fiber-dominated  properties.  Leaving  the  fiber  dominated 
properties  unchanged,  we  specifically  seek  to  develop  an  averaging  procedure  for  the  less 
important  three  matrix-dominated  properties.  The  procedure  is  intended  to  simplify  the 
constitutive  form,  but  it  must  not  significantly  alter  basic  physical  behavior.  Needless  to  cay, 
the  matrix  dominated  properties  averaging  procedure  will  be  of  no  interest  or  use  if  it  does  not 
permit  the  development  of  a  3-D  lamination  theory  and  the  subsequent  detailed  evaluation 
there  of.  Furthermore,  it  remains  possible  tha‘  the  entire  procedure  could  be  "exact"  in  special 
cases,  as  will  turn  out  to  be  true. 


2.  THE  FORMAL  METHOD 

The  macroscopic  elastic  properties  for  the  aligred  fiber  reinforced  medium  are  taken  to  be 
those  of  transversely  isotropic  symmetry  with 
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where  a*  =  C-fc  and  with 


and 


°1  -  °11*  °2  -  °22*  °3  -  °33>  °4  “  °23>  °5  “  °31>  °6  “  °12 
^1  =  ^n*  ^2  =  ®11'  ^3  =  ^33’  ^4  =  ^23’  %  =  ^31*  ^6  =  ^12  • 


The  five  independent  properties  in  (1),  with  axis  1  in  the  fiber  direction,  can  be  written  in  terms 
of  the  usual  measurable  properties,  E1V  v12,  £22,  n12.  and  ^3  through: 


Cn  =  Eu  +  4v12k23  , 
C12  =  2k23v12  , 


^22  =  M-23  +  K23  > 
^23  *  “^23  +  K23  • 

^66  =  f*12  • 


(2) 


Here 


4  [  1  v^  E22/E11]  E22/h23  , 

with  Eu  and  Ezz  being  uni-axial  moduli,  v12,  the  axial  Poisson's  ratio,  and  |j.12  and  p23  the 
shear  moduli. 

The  grouping  of  five  independent  properties,  Eu,  v12,  Ezz,  n12,  and  contains  the  three 
matrix-dominated  moduli,  Ezz,  412,  and  which  were  discussed  previously.  To  permit 
developing  a  formalism,  two  interrelations  are  taken  between  the  three  matrix  dominated 
properties.  After  obtaining  the  relevant  forms,  consideration  will  be  returned  to  material 
systems  for  which  the  two  interrelationships  do  not  apply,  which  is  the  case  of  primary 
interest. 

Begin  by  taking  the  two  shear  moduli,  p12  and  pjg,  as  being  equal,  i.e., 

M-1 2  =  ■  W 


3 


The  second  relationship  between  the  matrix-dominated  moduli  is  taken  as 

( 1  ^12)  E 22 


M^23  * 


2(1  -  vf2  E22/£11) 


(5) 


(This  form  actually  results  from  a  corresponding  relationship  between  the  C,y  moduli,  namely, 
C12  =  C23.)  Examining  modulus  Cn  in  light  of  the  two  forms  (4)  and  (5),  Cu  can  be 
expressed  as 

C11  *  EU  -  E  +  6,1  .  (6) 

where 


and 


E  -  „(1  ~Vl2)£z2  -  2(1  +v12)|i12  =  2(1  +  v12) |x23 
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With  ,  so  defined,  for  the  other  C„’s  take 


Aj  =  Cljt  when  /,  j  *  1 ,  1  . 


Then  the  full  vector  of  d’s  subject  to  (4)  and  (5)  can  be  written  as 
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The  form  (10)  rigorously  applies  only  when  (4)  and  (5)  are  satisfied,  but  it  does  provide  a 
guiding  formalism  in  the  more  general  case  of  interest  here. 


At  this  point,  consideration  is  returned  to  fiber  reinforced  systems  which  do  not  satisfy  the 
two  interrelationships,  (4)  and  (5),  and  accordingly  accommodate  five  independent  properties. 
For  such  systems,  it  is  our  intention  to  develop  an  averaging  procedure  for  E22,  p.12,  and 
which  depends  explicitly  on  all  three  properties.  Of  course  these  moduli  cannot  be  averaged 
directly  since  E22  is  a  uni-axial  modulus,  while  the  other  two  are  shear  moduli.  From  (10), 
however,  a  format  is  evident  by  which  E22  can  be  compared  with  p12  and  Define 
generalized  moduli  by: 

gj  0  ~  V12)  E 22 

1  "  2(1- v122E22/Eu)  ’ 

^2  *  M-12  ’ 

S*3  =  H23-  (11) 


An  averaged,  generalized  modulus,  S£,  could  be  established  from  (11)  by  arbitrarily  writing, 
S£  =  (s£,  +  S£2  +  S£a)/3.  However,  a  better  and  more  physically  meaningfully  averaging 
procedure  now  will  be  developed. 

To  establish  a  base-line  accessible  case  for  physical  interpretation,  examine  the  special 
form  of  transverse  isotropy  which  results  when  the  two  Poisson’s  ratios  are  taken  as 
vanishing,  i.e., 


v12  =  0  and  v23  =  0  . 


(12) 


From  (2)  it  immediately  follows  for  v12  =  0  that  C12  =  0.  The  identity  E22  =  2(1  +  v23)p23 
combined  with  (12)  produces  E22  =  2^3.  Using  this  latter  result  in  (3)  and  then  (2)  results  in 
C23  =  0.  Finally,  using  E22  =  2^3  again  in  (3)  and  (2)  gives  C22  =  E22,  so  that  (12)  with  (2) 
yields  Cn  =  E1V  Combining  all  these  results  into  (1)  gives  it  the  diagonal  form 
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E„  0  0  0  0  0 

e22  0  0  0  0 

E22  0  0  0 

H*3  0  0 

^2  0 
^2 

Thus,  (13),  with  E22  =  2^3,  is  the  form  assumed  by  transverse  isotropy  for  vanishing 
Poisson’s  ratios,  (12),  with  no  other  restrictions. 

Now  form  the  strain  energy  density,  U,  associated  with  (13).  This  gives 

2.U  -  (  E^  -  Ej2)  Eli  +  E22Cii  +  E22£22  +  £^£33  +  2|i12C12  +  2^23 e 23  +  2(112E31  .  (14) 

The  term  (En  -  E22)z2u  is  the  fiber-dominated  effect,  which  obviously  must  vanish  as  the  fiber 

reinforcement  contribution  goes  to  zero  in  an  isotropic  matrix  phase.  For  the  remaining 
matrix-dominated  terms  in  (14)  let 


U  -  U, 


Matrix -Dominated 


£?22  2  2  2  2  2  2 
-  — pj-  (*n  +  +  E33)  +  H23eZ3  +  p12(c12  +  e31)  . 


In  the  spirit  of  considering  fully  3-D  deformation  conditions  take  all  strain  components  as 
Ey  =  O  (5),  and  using  this  in  (15)  gives 

(v  1  \ 

0=O  3^i  +  2Hl2  +  n23  52  .  (16) 

U  2  -I  ) 

It  is  seen  that  the  three  matrix-dominated  moduli  contribute  to  the  strain  energy  (16)  in  the 
proportions 
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m  2  •  1^23 


as 


3:2:1. 


(17) 


Returning  to  the  generalized  moduli,  (11),  an  acceptable  generalized  averaging  procedure 
must  recover  the  result  just  found.  Specifically,  for  the  generalized  moduli  in  (1 1)  (specialized 
to  v12  =  v23  =  0)  to  recover  the  results  (16)  and  (17),  it  is  necessary  that  the  moduli  be 
averaged  with  the  weighting  factors  shown  as  follows: 

:  S£2  :  as  3:2:1.  (18) 

Thus,  this  generalized  averaging  procedure  gives  the  generalized  shear  modulus,  p,  as 


1 
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3(1  ^12)  ^22 

2(1  -v52eM/£„) 


+  2m2  +  ^23 


(19) 


Finally,  then  the  coefficients  formalized  by  (10)  are  now  taken  as 


1  - 
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(20) 
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where  |x  is  now  given  by  (19).  Relation  (20)  also  applies  in  the  "exact"  case  when  (4)  and  (5) 
are  identically  satisfied. 

Next,  write  the  £y’s  in  (20)  in  matrix  form  and  compare  them  with  the  appropriate  isotropic 
material  form  in  Green  and  Zema  (1963).  It  is  quickly  seen  that  the  forms  in  (20)  are  those 
for  isotropic  behavior  with  isotropic  shear  modulus  given  by  |i.  and  with  isotropic  X  given  by 


X  = 


2ViZ 

1  -  2v12 


V-  ■ 


(21) 


The  relation  between  £(j  and  Cy  is  given  by  (6)  and  (9).  The  form  is  isotropic,  as  just 
established.  Thus,  using  (9),  the  Ctj  form  is  isotropic  to  within  the  presence  of  the  (£n  -  £) 
part  of  C1r  It  directly  follows  that  the  stress  constitutive  relation,  corresponding  to  the  Cy 
form,  is  given  by 


°ij  =  +  2^j  +  (^n-  Wifi  i  .  (22) 

where  X  is  given  by  (21),  \i  by  (19)  and  £  is  given  by 


£  =  2(1  +  v12)n  •  (23) 

The  form  (22)  was  previously  identified  by  Christensen  (1988);  however,  in  that  work  the 
associated  interrelationships  (4)  and  (5)  were  taken  as  restrictions,  rather  than  being 
developed  into  the  generalized  averaging  procedure  given  here  and  culminating  in  expression 
(19).  Also  in  the  work  of  Christensen  (1988),  emphasis  was  given  to  using  (22)  to  develop  an 
associated  failure  criterion.  Attention  here  is  given  to  using  the  form  (22)  to  identify  a  specific 
3-D  lamination  procedure  and  to  evaluate  it. 

The  lamina  level  constitutive  formulation  is  now  complete  with  the  forms  (19)  and  (21-23). 

It  should  be  emphasized  that  the  lamina  constitutive  equation  (22)  does  not  imply  the 
interrelations  (4)  and  (5),  although  it  only  is  "exact"  when  they  are  satisfied.  In  the  most 
general  case,  all  five  transversely  isotropic  properties  enter  the  constitutive  form  (22).  The 
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two  fiber-dominated  properties  enter  the  constitutive  forms  directly,  while  all  three  matrix- 
dominated  properties  enter  through  the  generalized  averaging  procedure  result  (19).  Using 
the  form  (22)  to  specifying  the  lamina-level  behavior,  the  corresponding  laminated  medium 
characterization  will  be  now  found.  This  is  the  so  called  3-D  lamination  theory  or  procedure. 

Express  constitutive  relation  (22)  in  coordinates  other  than  in  the  fiber  direction.  Take  the 
angle  of  rotation  of  the  coordinate  system  as  9  from  the  fiber  direction,  rotated  about  axis  x3. 
With  3jj  being  the  direction  cosines  matrix,  the  constitutive  relation  (22)  takes  the  rotated  form 

Ojj  =  +  2p£jj  +  (EU-  E)auana^a^  ,  (24) 

where 

[a,J  =  [cos0,  -sine,  0] ,  (25) 

with  (25)  being  the  first  row  of  a^.  Now  take  N  similar  lamina  in  bonded  contact,  each  with  its 
individual  direction,  en.  The  3-D  constitutive  form  (24)  for  each  lamina  is  directly  combined  to 
give  the  3-D  laminated  medium  constitutive  form 


=  X.i 


kk  &ij  + 


2  HV 


n*  1 


N 


_(n)  _(«)  _("> 

m,  m.  mk  m,  z 


kl  ’ 


(26) 


wherein  from  (25) 


=  [cos6n,  -sine„,  0]  .  (27) 

Relation  (26)  along  with  (19),  (21),  (23),  and  (27)  constitute  the  3-D  constitutive  theory  for  fiber 
composite  laminated  media  assembled  from  a  single  lamina  type,  and  involves  all  five 
independent  properties. 

It  is  readily  verified  from  (26)  that  the  interlaminar  stress  components  o33,  o32,  and  o31  are 
continuous  across  the  interfaces  between  lamina  because  the  last  term  in  (26)  vanishes  for 
these  three  stress  components.  By  definition  the  displacements  are  continuous  across  the 
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interfaces  since  the  entire  constitutive  form  (26)  is  expressed  in  terms  of  the  single  strain 
tensor,  e^.  Thus,  the  laminated  medium  constitutive  form  (26),  as  derived  from  the  lamina 
constitutive  form  (22),  is  an  exact  3-D  result.  No  approximations  are  involved  in  satisfying 
interface  conditions.  It  is  implicit  in  obtaining  (26)  from  (22)  that  the  strain  gradients  are  small 
over  the  number  of  lamina  of  interest,  N.  This  condition  will  be  generalized  (relaxed)  in  the 
section  on  utilization.  Although  form  (26)  takes  all  lamina  to  be  of  equal  thickness,  it  is  easily 
extended  by  including  thickness  related  weighting  functions  in  it. 

The  constitutive  form  (26)  for  the  laminated  medium  is  extraordinarily  compact  and  easy  to 
use.  If  the  two  interrelations  (4)  and  (5)  are  satisfied  (to  within  experimental  accuracy)  by  a 
given  set  of  five  transversely  isotropic  properties,  then  the  entire  procedure  up  through  and 
including  the  final  constitutive  result,  (26),  is  "exact"  within  the  theory  of  linear  elasticity.  In 
this  case,  the  results  are  identical  to  the  Cj,  matrix  in  Pagano  (1974).  If  the  two  interrelations 
(4)  and  (5)  are  not  satisfied,  then  the  result,  (26),  is  an  approximation  based  upon  the 
generalized  averaging  procedure  for  the  three  matrix-dominated  properties.  Under  this 
condition  it  will  be  necessary  to  carefully  evaluate  the  nature  of  the  approximation  as  will  be 
done  in  the  next  two  sections. 


3.  A  TEST  CASE 

Judging  the  total  accuracy  of  the  proposed  constitutive  relationship  against  the  exact  form 
necessitates  comparing  more  than  just  the  resulting  effective  lamina  moduli.  Although  direct 
comparisons  between  individual  moduli  show  the  differences  expected  from  "uni-axial"  load 
states,  they  do  not  accurately  capture  the  coupling  between  components  which  occurs 
naturally  under  a  general  3-D  load  state.  This  section  formulates  a  specific  BV  problem  which 
first  allows  direct  comparisons  of  various  field  quantities  obtained  by  using  both  the  exact  and 
proposed  constitutive  relationships,  and  second  shows  the  finite  regimes  where  classical  long- 
wavelength  (LW)  theory  (essentially  plane  stress  or  classical  flat  lamination  theories)  is 
appropriate.  To  aid  the  reader  in  distinguishing  variables  associated  with  the  constitutive 
formulation  from  those  associated  with  the  BV  problem,  direction  notation  and  an  X-Y-Z 
coordinate  frame  are  employed  in  defining  and  solving  the  test  problem. 
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In  adopting  a  test  case  for  the  3-D  lamination  theory,  it  is  necessary  to  select  a  particular 
lay-up  pattern.  There  are  two  limiting  cases  for  lay-up  patterns.  At  one  extreme  is  the 
degenerate  case  of  all  lamina  being  aligned,  thus  the  laminate  retains  aligned-axis  character. 
At  the  other  extreme  is  a  quasi-isotropic  lay-up  involving  the  most  dispersed  pattern  of 
directions.  The  aligned  laminate  is  a  trivial  application  of  lamination  theory,  since  there  are  no 
true  interfaces.  The  quasi-isotropic  lay-up  is  in  fact,  the  most  extreme  case,  and  provides  the 
most  severe  test  of  a  lamination  theory.  It  requires  the  highest  degree  of  interactive  coupling 
between  all  lamina  properties  to  produce  laminate  behavior. 

3.1  The  Boundary  Value  Problem.  The  problem  considered  herein  is  a  "homogenized" 
solid  lying  parallel  to  the  X  -  Z  plane.  The  laminated  solid  consists  of  many  identical 
transversely-isotropic  laminae  oriented  in  a  repetitive  [0°,  ±60°]  order  and  is  considered  to  be 
quasi-isotropic  in  the  X  -  Z  plane.  A  sufficient  number  of  laminae  exist  through  the  thickness, 
in  the  Y  direction,  such  that  the  macroscopic  composite  properties  are  well  represented  by 
those  of  a  homogenized  single  [0°,  ±60°]  sublaminate. 

Without  loss  of  generality,  the  coordinate  positions  X,  Y,  and  Z  are  non-dimensionalized  to 
the  coordinates  x,  y,  and  z  by  half  the  thickness  of  the  solid,  hJ2,  ( i.e x  =  2X/h,  y  =  2Y/h,  and 
z  =  2Z/h)  so  that  the  upper  and  lower  surfaces  lie  in  the  planes  y  =  1  and  y  =  -1 ,  respectively. 
Displacements  in  the  z  direction,  as  well  as  all  derivatives  taken  with  respect  to  z,  are  taken  to 
vanish  identically  (yielding  plane  strain  conditions).  Imposed  on  the  lower  surface,  y  =  -1,  are 
traction-free  boundary  conditions;  while  on  the  upper  surface,  y  =  1 ,  a  sinusoidal  normal 
traction  is  applied  which  has  the  form 


tn  =  -P  sin(cu)  . 


(28) 


Here 


(29) 


P,  the  load  magnitude,  has  units  of  stress,  2 L  is  the  loading  wavelength,  and  the  convention 
that  positive  normal  tractions  produce  normal  tensile  stresses  is  employed.  The  solid  extends 
"infinitely"  in  the  x  direction,  but  only  the  region  from  x  -  0  to  x  =  2L/h  need  be  examined,  due 
to  periodicity.  Furthermore,  all  body  forces  are  identically  zero.  This  problem  has  been 
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chosen  because  by  varying  a,  essentially  the  ratio  of  h  to  L,  many  different  loading  conditions 
can  be  achieved.  For  values  of  h/L  «  1 ,  classical  LW  conditions  result,  and,  as  h/L 
increases,  fully  3-D  load  states  develop. 

3.2  Material  Formulation.  Using  the  homogenizing  procedure  used  by  Sun  and  Li  (1988), 
the  effective  sublaminate  stiffnesses  are  assembled.  Using  the  appropriate  tensor 
transformations,  the  ±60°  material  stiffness  tensor  components  are  expressed  in  the  0°  lamina 
material  coordinate  frame.  (For  transversely  isotropic  lamina  whose  plane  of  isotropy  lies 
perpendicular  to  the  lamina-  or  ply-  plane,  Christensen  (1988)  gives  the  explicit  stiffness 
tensor  components  in  any  cartesian  coordinate  frame  where  the  rotated  inplane  axes  remain 
parallel  to  the  original  ply-plane.)  Unit  macro  strains  are  individually  imposed  along  the 
boundaries  of  a  combined  (unit-square)  ±60°  lamina  pair.  After  enforcing  continuity  of 
interlaminar  displacements  and  tractions,  the  resulting  net  forces  yield  the  effective  ±60° 
stiffness  components  when  properly  normalized.  Next,  the  0°  and  effective  ±60°  laminae  are 
weighted  appropriately  and  assembled  using  the  same  method  to  produce  the  effective 
properties  of  the  entire  sublaminate. 


The  quasi-isotropic  lay-up  produces  a  set  of  transversely  isotropic  properties  in  three 
dimensions.  With  x3  being  the  axis  of  symmetry  for  the  laminate,  the  non-zero  effective 
stiffness  properties,  cljt  for  the  [0°,  ±60°]  quasi-isotropic  laminate  are  given  by  these  newly 
derived  (exact)  closed  form  expressions: 


C11  =  C22  =  •g  (^11  +  ^22)  +  ^12  +  ^66  QQ~  ^  ~  ^23^  ’ 

i <C„  *  c!2)  ♦  1  C„-  C„ .  ’  <C„ -  . 

^13  *  ^23  “  ^  (^12  *  ^23)  ' 

C33  *  ^22  ' 


C*.  -  Css  - 


^  C(66  (  ^22  “  ^23  ) 


66  \  22 
2  Cgj  +  C22 


r 23 


■ 2  (  cn  ciz)  • 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 
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Here  the  upper  case  C(j  refers  to  the  individual  lamina  components  given  in  the  lamina 
coordinate  frame;  i.e.,  axis  1  in  the  fiber  direction.  Substituting  (2)  and  (3)  into  (30-35)  allows 
the  quasi-isotropic  stiffnesses  to  be  expressed  in  terms  of  the  usual  five  properties,  specifically 
£|1>  ^22'  v12’  M-12'  ^nd  ^23' 

3.3  Solution  Formulation.  The  general  solution  to  the  test-case  BV  problem  has  been 
fully  derived,  following  the  method  employed  by  Timoshenko  and  Goodier  (1970)  to  solve  a 
similar  isotropic  problem.  First,  the  equilibrium  equations,  infinitesimal  strain-displacement 
relationships,  and  compatibility  equations  are  reduced  by  incorporating  the  plane  strain 
assumption.  Next,  the  individual  stress  component  forms  are  found  by  using  the  reduced 
equilibrium  equations  and  by  assuming  a  separable  solution  where  a  Fourier  series  represents 
the  x  dependence  and  some  function  of  F(y)  represents  the  y  dependence.  The  strain 
component  forms,  determined  by  substituting  the  stress  forms  into  the  constitutive  relationship, 
along  with  the  sole  remaining  compatibility  equation  yield  the  governing  differential  equations 
for  F(y).  Solving  the  linear  homogeneous  differential  equation  and  matching  the  original 
boundary  conditions  uniquely  determines  F{y)  which,  when  substituted  back  into  the  original 
stress  and  strain  forms,  yields  the  complete  closed-form  BV  problem  solution. 

The  effective  laminate  properties  and  BV  problem,  in  conjunction  with  actual  material 
properties,  allow  an  objective  and  realistic  evaluation  of  the  proposed  constitutive  relationship 
in  a  fully  3-D  load  state.  In  the  following  section,  laminate  properties  and  solutions  to  the  BV 
problem  will  be  compared  for  two  specific  composite  systems,  and  the  effectiveness  of  the 
proposed  lamination  procedure  relationships  will  be  demonstrated. 


4.  EVALUATION 

4.1  Properties.  To  evaluate  the  proposed  constitutive  relationship,  two  sets  of  lamina 
properties  are  chosen  which  approximately  satisfy  the  restrictions  of  transverse  isotropy.  A 
graphite/epoxy  (Gr/Ep)  composite  [AS4/3501  -6]  is  chosen  for  its  extremely  anisotropic  (axial  to 
transverse)  moduli,  along  with  a  glass/epoxy  (Gl/Ep)  composite,  chosen  as  representative  of  a 
lesser  fiber-dominated  system.  Table  1  lists  the  elastic  lamina  constants  for  both  systems, 


13 


Table  1.  Lamina  Properties  Used  in  Evaluation 


Gr/Ep  t 

Gl/Ep  tt 

Eu  GPa 

144 

53.7 

^22  =  ^33  GPa 

9.65 

17.9 

V12 

0.31 

0.25 

V23 

0.52 

0.25 

|i12  GPa 

5.24 

8.96 

M-23  GPa 

3.17 

7.17 

t  Data  from  Kim,  Abrams,  and  Knight,  1988. 
tt  Data  from  Sun  and  Li,  1 988. 


and  to  satisfy  the  restrictions  of  transverse  isotropy,  namely  1^3  =  E22/ 2(1  +  v23),  the  actual 
moduli  are  adjusted  slightly  from  the  reported  values. 

Table  2  contains  the  effective  quasi-isotropic  laminate  properties  calculated  via  (30-35) 
using  the  proposed  forms  from  Section  2,  namely  the  C()  properties  corresponding  to  the 
lamina  constitutive  equation  (22),  and,  for  comparison,  using  the  exact  C,,  properties  from  (2) 
and  (3).  [The  proposed  Ctl  properties  are  explicitly  given  by  (6),  (9),  (19-21)  and  (23).]  The 
ratios  of  these  results  are  shown  in  Table  2  as  well.  Generally,  less  discrepancy  exists 
between  the  inplane  moduli  (1  and  2  directions)  than  the  out-of-plane  components.  Also 
tabulated  in  Table  2  are  the  quantities  En/(  1  -vf2)  and  En/(  1  -  v12)  which  represent  uni-axial 

and  bi-axial  inplane  extensional  stiffnesses,  respectively,  associated  with  classical  flat 
lamination  theory.  For  the  two  composite  systems  examined,  the  exact  and  proposed  quasi¬ 
isotropic  moduli  and  inplane  extensional  stiffnesses  differ  typically  by  less  than  =  10%. 

4.2  Test  BV  Problem.  The  previously  formulated  BV  problem  is  now  explored,  utilizing 
the  Gr/Ep  and  Gl/Ep  laminate  properties  listed  in  Table  2,  to  quantify  results  from  a  3-D 
problem.  By  judiciously  selecting  several  variables  at  two  locations,  a  meaningful  and 
representative  comparison  is  possible  without  having  to  examine  the  entire  field.  The 
variables  chosen  for  comparison  are:  first  the  flexural  stress  and  strain,  oxx(y  =  1,  x  =  Uh) 
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.  able  2.  Effective  Quasi-Isotropic  Laminate  Elastic  Constants  and  Inplane  Extensional 
Stiffnesses  Calculated  Using  the  Exact  Transversely  Isotropic  and  Proposed 
Lamina  Constitutive  Relationships  for  Gr/Ep  and  Gl/Ep  Systems 


Gr/Ep 

i 

Gl/Ep  j 

Proposed 

Exact 

ProDOsed 

Exact 

Proposed 

Exact 

Proposed 

Exact 

En  -  E22  GPa 

55.0 

55.5 

0.99 

30.7 

31.1 

0.99 

^33  GPa 

13.4 

12.5 

1.07 

20.4 

18.3 

1.11 

V12 

0.329 

0.309 

1.06 

0.283 

0.245 

1.16 

V13  =  v23 

0.302 

0.348 

0.87 

0.239 

0.225 

1.06 

p12  GPa 

20.7 

21.2 

0.98 

11.9 

12.5 

0.95 

^13  =  1^23  GPa 

3.95 

3.95 

1.00 

7.61 

8.00 

0.95 

£ii/(1-Vi22)  GPa 

61.7 

61.3 

1.01 

33.4 

33.1 

1.01 

Ei^l-v,,)  GPa 

82.0 

80.0 

1.02 

42.8 

41.1 

1.04 

and  e^y  =  0,  x  =  L/h);  second  the  upper-surface  normal  displacement.  uy(y  =  1,  x  =  L/h);  and 
third  the  mid-plane  shear  stress  and  shear  strain.  oxy(y  =  0;  x  =  L/h)  and  cxy(y  =  0,  x  =  L/h). 

At  these  locations,  these  variables  are  the  dominant  components  present  as  well  as  the 
minimum  or  maximum  values  achieved  by  these  variables  (almost)  anywhere  in  the  body. 

The  normal  displacement  probably  best  captures  the  3-D  aspects  in  that  it  reflects  an 
integrated  value,  where  as  the  other  terms  are  more  indicative  of  the  material  properties 
directly  associated  with  that  particular  component. 

Table  3  lists  the  normalized  flexural  stresses,  as  a  function  of  the  thickness-to-length  ratio, 
h/L,  evaluated  using  both  the  proposed  and  exact  constitutive  relationships  for  both  composite 
systems.  The  stresses  have  been  normalized  by  the  asymptotic  solution  obtained  as  h/L-*  0, 
i.e.,  essentially  LW  or  classical  plate  theory,  using  the  exact  transversely-isotropic  lamina 
relationship.  Reasonable  agreement  exists  between  the  two  solutions  from  h/L  =  0  to  an 
astonishing  value  of  h/L=  10.  Note  how  rapidly  the  BV  solution  deviates  from  classical  LW 
theory  as  h/L  increases  beyond  approximately  0.10.  Figures  1  and  2  show  the  mid-plane 
shear  strain  and  flexural  strain,  respectively,  normalized  by  their  LW  solution  for  both 
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Table  3.  Normalized  Flexural  Stresses  for  Gr/Ep  and  Gl/Ep  Systems 
Verses  Thickness  to  Length  Ratio  h/L 


h/L 

Gr/Ep 

Gl/Ep  j 

Proposed 

Exact 

Proposed 

Exact 

0.01 

1.000 

1.000 

1.000 

1.000 

0.05 

1.006 

1.006 

1.001 

1.001 

0.10 

1.023 

1.023 

1.006 

1.006 

0.50 

1.663 

1.638 

1.221 

1.216 

1.00 

3.824 

3.965 

2.319 

2.406 

5.00 

86.18 

90.39 

51.60 

54.78 

10.0 

343.1 

359.8 

205.4 

218.1 

Note:  Normalized  flexural  sti  jsses  are  defined  by  axx  (y  -  1,  x  -  L/h)/axx.LW(y  -  1,  x  -  Uh), 
where  the  normalizing  value  axx.LW,  represent  the  asymptotic  solution  obtained  as  h/L  -»  0; 
i.e.,  LW  theory. 


composite  systems  as  a  function  of  h/L  The  agree  lent  between  the  results  shown  in 
Figures  1  and  2  is  typical  of  all  the  variables  considered.  Examining  the  figures  reveals  that 
the  components  calculated  from  the  proposed  and  exact  constitutive  relationship  differ  only 
slightly  even  as  -*  0  as  h/L  ->  °°.  Unfortunately,  plotted  at  this  scale  it  is  not  evident  how 
drastically  the  BV  and  LW  solutions  drift  apart,  e.g.,  exx  >  2exx.LW  for  h/L  >  0.75  and 
Ejy  <  OJSe^ ,LW  for  h/L  >  1 .0.  Nonetheless,  the  proposed  lamina  relationship  captures 
remarkably  well  all  characteristics  of  the  BV  problem,  neglecting  small  differences  in  numerical 
values. 


To  quantify  precisely  the  differences  arising  from  the  proposed  relationship,  the  five  field 
quantities  have  been  normalized  by  their  corresponding  exact  valjes  and  plotted  verses  h/L  in 
Figures  3  and  4  for  the  Gr/Ep  and  Gl/Ep  systems,  respectively.  Recognize  that  as  h/L 
oxy  0  and  -»  0,  and  thus  their  normalized  ratios  lor  increasing  h/L  represent  differences 
in  decreasingly  smaller  numbers.  Ignoring  the  regions  where  a ^  and  exy  approach  zero,  the 
errors  in  the  five  normalized  field  quantities  for  h/L  <  10  are  less  than  1 1%  in  the  Gl/Ep 
system  and  7%  in  the  fiber-dominated  Gr/Ep  system. 
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Overall,  the  BV  problem  results  invariably  reflect  the  differences  in  effective  lamina 
properties  between  the  proposed  and  exact  constitutive  idealization.  Keeping  in  mind  that  tne 
BV  problem  is  intrinsically  stress  controlled,  underestimated  properties,  like  Gl/Ep  413,  produce 
overestimated  strain  values  and  vice  versa. 

This  test  problem  is  ideally  suited  to  discriminate  thin  laminate  conditions  (LW  theory)  from 
thick  laminate  conditions.  The  ratio  of  laminate  thickness  to  the  loading  wavelength  (h/L) 
directly  gives  the  parameter  of  variation.  As  h/L  -»  0,  long  wavelength,  thin  laminate 
conditions  are  recovered  as  a  limiting  case.  As  h/L  becomes  larger,  more  complex  physical 
effects  enter  the  BV  problem.  Somewhere  in  the  range  0.1  <  h/L  <  1  a  transition  from  thin 
laminate  conditions  to  thick  laminate  conditions  occurs,  e.g.,  see  Figures  1  and  2.  Certainly, 
in  the  neighborhood  of  h/L  =  1  fully  3-D  deformation  conditions  exist  in  the  laminate.  The 
lamination  theory  developed  herein  properly  models  all  major  effects,  and  it  even  remains 
effective  up  through  h/L  =  10. 

4.3  Constitutive  Error  Relative  to  Experimental  Uncertainties.  As  stated  in  Section  2,  this 
work  aims  to  formulate  a  practical  and  workable  3-D  lamination  theory  which  retains  algebraic 
simplicity  while  providing  "reasonable  and  useful"  composite  representation.  In  actual 
application  both  constitutive  approximations  and  material  property  variations,  arising  from 
experimental  uncertainties,  contribute  to  overall  analysis  inaccuracies.  Thus,  comparing  the 
absolute  mathematical  error  from  the  constitutive  approximations  with  that  from  the  material 
properties  appears  appropriate.  This  section  compares  the  magnitude  of  experimental 
uncertainties  to  the  resulting  constitutive  error  in  quantifying  the  lamina  properties. 

Large  scatter  in  measured  iamina  properties  is  common,  especially  with  all  the  difficulties 
encountered  in  obtaining  them.  For  example,  Kim,  Abrams,  and  Knight  (1988)  report 
coefficients  of  variations  in  measured  uni-directional  Gr/Ep  lamina  moduli  between  2.0%  and 
12.2%  even  though  3  to  10  specimens  were  used  in  each  test.  Clearly  the  lamina  constants 
most  difficult  to  experimentally  determine  and  most  prone  to  uncertainties  are  the  three  matrix- 
dominated  moduli,  namely,  E22,  p12,  and 

The  differences  between  the  proposed  and  exact  lamina  relationships  stem  solely  from 
assumptions  made  regarding  the  matrix-dominated  properties.  Since  variations  in  fiber- 
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dominated  properties  produce  (nearly)  identical  results  in  either  constitutive  relationship, 
attention  is  focused  only  upon  uncertainties  associated  with  the  matrix-dominated  properties. 
Plotted  in  Figures  5  and  6  are  the  normalized  BV  problem  quantities  for  the  Gr/Ep  and  Gl/Ep 
systems,  respectively,  calculated  by  decreasing  the  three  matrix-dominated  lamina  properties 
by  10%  (an  amount  typical  of  their  actual  variation)  from  their  base  values  given  in  Table  1. 
As  before,  the  BV  problem  values  are  normalized  by  their  corresponding  exact  quantities. 
The  variations  in  field  quantities  and,  thus,  lamina  properties  caused  by  uncertainties  in  the 
three  matrix-dominated  properties  are  the  same  order  of  magnitude  as  those  introduced  by 
the  proposed  constitutive  relationship;  i.e.,  compare  Figures  4  and  5  with  Figures  2  and  3. 
Therefore,  we  conclude  that  use  of  the  proposed  method  in  actual  engineering  analyses 
should  produce  results  with  uncertainties  which,  for  all  practical  purposes,  cannot  be 
realistically  differentiated  from  inherent  experimental  error  in  the  matrix-dominated  lamina 
properties. 


5.  UTILIZATION 

This  section  explores  how  the  proposed  constitutive  relationship  can  be  beneficially  utilized 
in  analyzing  both  thin  and  thick  composite  structures. 

5.1  Classical  Lamination  Theory.  The  proposed  lamina  relationship  can  be  incorporated 
into  classical  lamination  theory.  Using  the  previous  averaging  approach  on  the  four  constants 
required  by  classical  lamination  theory,  the  generalized  shear  modulus  becomes 


3(1  -v12)  E22 
2(1 -W ’2E22/Eu) 


(36) 


Using  (6),  (9),  (20),  and  (36),  the  transformed  n-th  lamina  stiffness  matrix  Q defined  such 
that 
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is  directly  obtainable.  Because  of  the  lamina  relationship  form,  Q,"  is  decomposable  into  an 
effectively  isotropic  and  homogeneous  term,  0"hH,  and  a  reinforcement  term,  Q,"_n>  as 

oj-  o,;.H+o,;.R,  (38) 


where 

[<v„r  « 
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0  Op 


(39) 


and 


[<V„r  *  (Eu-E) 
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Since  does  not  enter  into  classical  lamination  theory,  the  form  of  |x  given  in  (36)  should  be 
used.  When  evaluating  the  extensional,  coupling,  and  bending  stiffness  matrices  for  a 
laminate  composed  of  a  single  lamina  type,  the  contribution  of  (39)  factors  through  the 
integration,  as  a  homogeneous  isotropic  material,  leaving  only  the  reinforcement  contribution 
to  be  calculated  on  a  ply  by  ply  basis. 

5.2  Finite  Element  Procedures  for  Thick  Composite  Structures.  Currently,  to  analyze  a 
thick  composite  structure  either  each  lamina  in  the  component  is  discretized  and  assigned  its 
own  material  constants  or  the  laminae  are  homogenized,  represented  by  a  single  or  several 
sets  of  effective  properties,  allowing  the  component  to  be  discretized  as  an  equivalent 
"homogeneous"  solid.  When  defining  each  lamina  as  its  own  continuum  layer,  obtaining  a 
satisfactory  solution  with  modest  to  'arge  through  thickness  gradients  generally  requires 
several  elements  through  each  lamina  thickness;  and  because  typical  composite  structures 
contain  tens  to  hundreds  of  plies  and  because  element  aspect  ratios  must  be  maintained,  the 
total  Droblem  size,  i.e.,  degrees  of  freedom  (DOF),  quickly  escalates.  In  structures  where 
adjacent  parts  mandate  non-linear,  presumably  iterative,  solution  procedures  and  contain 
thousands  to  hundreds  of  thousands  of  DOF,  implicit  FE  methods  become  economically 
unfeasible.  On  the  other  hand,  when  dynamic  responses  necessitate  explicit  FE  schemes, 
e.g.,  impact  loadings,  the  maximum  time  increment  permitted  is  limited  by  approximately  the 
minimum  time  required  for  a  wave  to  propagate  across  the  smallest  element  (Bathe  1982). 
With  such  small  elements,  an  unacceptable  large  number  of  time  steps  results  and  may 
render  the  solution  technique  impractical. 

To  circumvent  some  of  these  difficulties,  consider  incorporating  the  lamina  constitutive 
form  (24)  into  displacement-based  solid  elements,  using  conventional  FE  methodologies. 
Following  the  usual  procedures,  the  laminated  structure  to  be  analyzed  would  be  discretized 
as  a  homogeneous  solid  (presuming  it  is  manufactured  from  a  single  composite  system). 
However,  when  evaluating  the  elemental  stiffness  matrix  K)y,  Ci}  is  allowed  to  vary  with  position 
representing  the  differently  oriented  lamina  within  the  element.  Symbolically,  K)y  is  thus 
calculated  as 

K'i  =  iB>lCtm(x,y.z)Bmjdv  ,  (41) 
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where  v  designates  the  element  volume,  C,j[x,  y,z)  is  the  position  dependent  stiffness  matrix 
whose  form  is  given  by  (24),  and  is  the  strain-displacement  transformation  matrix.  Note 
that  (41)  homogenizes  the  laminate  material  within  each  element  to  the  same  kinematic  order 
as  the  element,  without  sacrificing  stacking  sequence  related  behavior.  This  procedure 
guarantees  continuity  of  interlaminar  displacements  and  tractions  while  ensuring  that  both  the 
global  solution  and  homoyenizaiicr.  method  converge,  in  the  FE  sense,  as  the  number  of 
elements  increases.  Simplifications  in  calculating  (41)  can  be  made  by  requiring  the  elements 
be  oriented  so  that  one  surface  is  parallel  to  the  plane  of  the  lamina.  Although  the  resulting 
strains  would  be  the  individual  lamina  strains,  some  post-processing  is  necessary  to  recover 
the  inplane  lamina  stresses.  In  general,  this  approach  should  make  analyzing  thick  laminates 
involve  nearly  the  same  degree  of  (analysis)  difficulties  as  encountered  in  analyzing  thick 
homogeneous  solids. 


6.  CONCLUSION 

A  3-D  constitutive  theory,  which  rigorously  enforces  continuity  of  interlaminar  tractions  and 
displacements,  is  developed  for  thick  laminated  media  and  is  evaluated  by  direct  comparisons 
with  exact  solutions.  The  lamina  constitutive  relationship,  almost  isotropic  in  mathematical 
form  but  arbitrarily  anisotropic  in  physical  content,  simplifies  lamination  theory  without 
sacrificing  physical  meaning  or  even  significant  accuracy.  The  accuracy  level  of  this  theory  is 
well  inside  the  range  of  experimental  uncertainty  contained  in  the  properties  themselves. 
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Figure  2.  Flexural  Strains.  Normalized  bv  the  Exact  Asymptotic  Solution,  vs.  Thickness  to 


Loading  Wavelength  Ratio  for  a  Quasi-Isotropic  Laminate. 
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Proposed/Exact 


Figure  4.  Normalized  Field  Quantities,  as  a  Function  of  the  Thickness  to  Loading  Wavelength 


Ratio.  Calculated  Using  the  Proposed  Constitutive  Relationship  and  Normalized  b 


Their  Corresponding  Exact  Values,  for  Gl/E 
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Figures.  Field-Quantities,  Normalized  by  the  Exact  Solution,  vs.  the  Thickness  to  1  parting 
Wavelength  Ratio.  Values  Obtained  Using  the  Exact  Constitutive  Relationship  but 
by  Decreasing  the  Three  Matrix-Dominated  Properties  bv  10%.  for  Gr/Ep. 


Figure  6.  Field  Quantities,  Normalized  by  the  Exact  Solution,  vs.  the  Thickness  to  Loading 
Wavelength  Ratio.  Values  Obtained  Using  the  Exact  Constitutive  Relationship  but 
by  Decreasing  the  Three  Matrix  Dominated  Properties  bv  1 0%,  for  Gl/En 
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